An improved algorithm for computing altitude dependent corrected geomagnetic (CGM) coordinates from geocentric coordinates (and, where it exists, the inverse) using spherical harmonics is described. The method uses a tenth order spherical harmonic fit to the direction cosines (a unit vector) in a suitably chosen intermediate, altitude adjusted coordinate system. The need for this auxiliary coordinate system is to avoid convergence problems associated with the discontinuity in the CGM latitude at the magnetic equator at non- zero altitude. Altitude dependence is obtained by computing the spherical harmonic fits to CGM (and inverse) at 0, 300 and 1200 km altitude, and using a quadratic fit to interpolate each coefficient. The new algorithm provides a good representation of the CGM compression around the South Atlantic Anomaly in addition to modeling the increasing discontinuity with altitude at the magnetic equator. Comparisons are provided with previous approaches. Accuracy limitations and consistency between the direct and inverse computations are also discussed.